Integration across different omics datasets
Andrej Blejec
28. 4. 2024
1 Introduction
[Needs more text]
In this document we present a work-flow for integration across different omics datasets.
[Note] This is not the final version of the document.
2 File management system
Files are managed under the pISA-tree framework. To take advantage of pISA-tree metadata information, we will use package pisar (For installation from GitHub see pisar).
library(pisar)
pisa <- pisa()
# Set output directory (needed for .Rmd -> HTML)
pisa$oroot <- file.path(.aroot,"output")
.oroot <- pisa$oroot%%{r initialize,echo=FALSE,include=FALSE,message=FALSE,warning=FALSE} %%options(width=70) %%#library(Hmisc) %%#library(amisc) %%#library(Biobase) %%#library(xlsReadWrite) %%#library(xtable) %%#library(MASS) %%catln <- function(...) cat(...,"\n") %%
Final document is in subdirectory reports:
basename(outputFile)## [1] ".pdf"2.1 Info from pISA structure
Input/data directory
.inroot## [1] "../../../input"Results directory
.oroot## [1] "../../../output"## project: _p_Omics## Investigation: _I_Omics## Study: _S_multiOmics## Assay: _A_multiOmics-integration-R3 Packages and functions
%For canonical correlation analysis we will use package , which is used in our projects. %{r message=FALSE,warning=FALSE} %library(CCA) %
A package with regularized CCA and multiomics DIABLO method is mixOmics. Package igraph is needed for network analysis.
library(mixOmics)
library(igraph)Package for complex heatmaps.
#BiocManager::install("ComplexHeatmap")
library(ComplexHeatmap)
# https://www.rdocumentation.org/packages/pheatmap/versions/1.0.12/topics/pheatmap
library(pheatmap)3.1 Functions
Some original and adapted functions can be found in the file that is silently processed here.
%% Additional functions
out <- ""
out <- paste(out,knit_child("005-Functions.Rmd", quiet=TRUE))4 Data
It is advisable to first read phenodata and featuredata, followed by actual data input. This enables smart selection of samples, based on the sample selection column with the assay name.
4.1 Phenodata
(pfn <- getMeta(.ameta
, "Phenodata"))## [1] "../../phenodata_20221001.txt"dir(file.path(.aroot,dirname(pfn)), pattern = basename(pfn))## [1] "phenodata_20221001.txt"phdata <- read.table(file.path(.aroot,pfn)
, header = TRUE
, sep = "\t"
, stringsAsFactors = FALSE
, row.names=1
)
dim(phdata)## [1] 32 15names(phdata)## [1] "SampleID" "Treatment"
## [3] "Harvest" "SamplingDay"
## [5] "DaysOfStressH" "PlantNo"
## [7] "Sample.type" "Date"
## [9] "Heat.Recovery.Days" "TreatmentxDatexPlant"
## [11] "TreatmentxSamplingDay" "TreatmentxSamplingDayxPlantNo"
## [13] "Transcriptomics" "Metabolomics"
## [15] "Hormonomics"Check and use the sample selection column, if present.
.aname## [1] "_A_multiOmics-integration-R"selectId <- substr(gsub("-",".",.aname),2,nchar(.aname))
selectId <- .vzorci
selectId## [1] NAif(selectId %in% names(phdata)) pdata <- phdata[!is.na(phdata[,selectId]),] else
pdata <- phdata
dim(pdata)## [1] 32 15Overview of selected samples:
table(pdata$Treatment, pdata$SamplingDay)##
## 1 7 8 14
## C 4 4 4 4
## H 4 4 4 4.treat <- unique(pdata$Treatment)
.days <- unique(pdata$SamplingDay)
.entry <- 0.5summary(pdata)## SampleID Treatment Harvest SamplingDay
## Length:32 Length:32 Min. :1.00 Min. : 1.0
## Class :character Class :character 1st Qu.:1.75 1st Qu.: 5.5
## Mode :character Mode :character Median :2.50 Median : 7.5
## Mean :2.50 Mean : 7.5
## 3rd Qu.:3.25 3rd Qu.: 9.5
## Max. :4.00 Max. :14.0
## DaysOfStressH PlantNo Sample.type
## Min. : 0.00 Min. : 7.00 Length:32
## 1st Qu.: 0.00 1st Qu.:10.75 Class :character
## Median : 0.50 Median :14.50 Mode :character
## Mean : 3.75 Mean :14.50
## 3rd Qu.: 7.25 3rd Qu.:18.25
## Max. :14.00 Max. :22.00
## Date Heat.Recovery.Days TreatmentxDatexPlant
## Length:32 Length:32 Length:32
## Class :character Class :character Class :character
## Mode :character Mode :character Mode :character
##
##
##
## TreatmentxSamplingDay TreatmentxSamplingDayxPlantNo Transcriptomics
## Length:32 Length:32 Min. :1
## Class :character Class :character 1st Qu.:1
## Mode :character Mode :character Median :1
## Mean :1
## 3rd Qu.:1
## Max. :1
## Metabolomics Hormonomics
## Min. :1 Min. :1
## 1st Qu.:1 1st Qu.:1
## Median :1 Median :1
## Mean :1 Mean :1
## 3rd Qu.:1 3rd Qu.:1
## Max. :1 Max. :1apply(pdata,2,function(x) summary(as.factor(x)))## $SampleID
## AD001 AD002 AD003 AD004 AD005 AD006 AD007 AD008 AD013 AD014 AD015
## 1 1 1 1 1 1 1 1 1 1 1
## AD016 AD017 AD018 AD019 AD020 AD025 AD026 AD027 AD028 AD037 AD038
## 1 1 1 1 1 1 1 1 1 1 1
## AD039 AD040 AD045 AD046 AD047 AD048 AD057 AD058 AD059 AD060
## 1 1 1 1 1 1 1 1 1 1
##
## $Treatment
## C H
## 16 16
##
## $Harvest
## 1 2 3 4
## 8 8 8 8
##
## $SamplingDay
## 1 7 8 14
## 8 8 8 8
##
## $DaysOfStressH
## 0 1 7 8 14
## 16 4 4 4 4
##
## $PlantNo
## 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
## 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
##
## $Sample.type
## adult leaf
## 32
##
## $Date
## 04/11/2020 10/11/2020 11/11/2020 17/11/2020
## 8 8 8 8
##
## $Heat.Recovery.Days
## 0_0 1_0 14_0 7_0 8_0
## 16 4 4 4 4
##
## $TreatmentxDatexPlant
## C_2020-11-04_10 C_2020-11-04_7 C_2020-11-04_8 C_2020-11-04_9
## 1 1 1 1
## C_2020-11-10_11 C_2020-11-10_12 C_2020-11-10_13 C_2020-11-10_14
## 1 1 1 1
## C_2020-11-11_15 C_2020-11-11_16 C_2020-11-11_17 C_2020-11-11_18
## 1 1 1 1
## C_2020-11-17_19 C_2020-11-17_20 C_2020-11-17_21 C_2020-11-17_22
## 1 1 1 1
## H_2020-11-04_10 H_2020-11-04_7 H_2020-11-04_8 H_2020-11-04_9
## 1 1 1 1
## H_2020-11-10_11 H_2020-11-10_12 H_2020-11-10_13 H_2020-11-10_14
## 1 1 1 1
## H_2020-11-11_15 H_2020-11-11_16 H_2020-11-11_17 H_2020-11-11_18
## 1 1 1 1
## H_2020-11-17_19 H_2020-11-17_20 H_2020-11-17_21 H_2020-11-17_22
## 1 1 1 1
##
## $TreatmentxSamplingDay
## C_1 C_14 C_7 C_8 H_1 H_14 H_7 H_8
## 4 4 4 4 4 4 4 4
##
## $TreatmentxSamplingDayxPlantNo
## C_1_10 C_1_7 C_1_8 C_1_9 C_14_19 C_14_20 C_14_21 C_14_22
## 1 1 1 1 1 1 1 1
## C_7_11 C_7_12 C_7_13 C_7_14 C_8_15 C_8_16 C_8_17 C_8_18
## 1 1 1 1 1 1 1 1
## H_1_10 H_1_7 H_1_8 H_1_9 H_14_19 H_14_20 H_14_21 H_14_22
## 1 1 1 1 1 1 1 1
## H_7_11 H_7_12 H_7_13 H_7_14 H_8_15 H_8_16 H_8_17 H_8_18
## 1 1 1 1 1 1 1 1
##
## $Transcriptomics
## 1
## 32
##
## $Metabolomics
## 1
## 32
##
## $Hormonomics
## 1
## 32%%## Featuredata %% %%{r } %%(ffn <- getMeta(.ameta,"Featuredata")) %% %% %%{r } %%if(ffn != "") %%fdata <- read.table(file.path(.iroot,ffn) %% , sep = "\t" %% , header = TRUE %% , na.strings = c("", "-") %% , stringsAsFactors = FALSE %% , row.names = 1 %% ) else fdata <- NULL %%head(fdata) %% %% %%First few columns, if present. %% %%{r } %%fdata[,1:2] %%
4.2 Process data files
For this project we aim to integrate several multi-omics datasets. We have data on hormonomics, metabolomics, and qPCR:
%%In the assay metadata file keys for data files start with keyword ‘Data’. %%Extract object names from those lines. %% %%{r} %%ind <- grep("^Data",.ameta[,1]) %%datanames <- gsub("^Data (.*):$","\\1",.ameta[ind,]$Key) %%datanames %% %%
%% %% %%Read data and put them in variables named indatasets: %%{r} %%cat("Created objects:\n\n") %% %%for (i in ind) { %%nm <- gsub("^Data (.*):$","\\1",.ameta[i,]$Key) %%x <- read.table(file.path(.aroot,.ameta[i,]$Value), header=TRUE, sep="\t") %%assign(nm,x) %%catln("\nDataset:", ls(pattern=nm)) %%str(get(nm)) %%} %%
hormonomics <- read.table(file.path(.aroot,"./input/data_hormonomics.txt"), header=TRUE, sep="\t")
metabolomics <- read.table(file.path(.aroot,"./input/data_metabolomics.txt"), header=TRUE, sep="\t")
qPCR <- read.table(file.path(.aroot,"./input/data_qPCR.txt"), header=TRUE, sep="\t")For future use, we need the text names of dataset objects:
datanames <- c("hormonomics", "metabolomics", "qPCR")Datasets for DIABLO need to be collected in a list, with rows corresponding to the same samples. The order of samples from shrinked phenodata will be enforced.
The first component of the list will be a grouping variable, indicating the conditions. We will create reasonable names for groups.
sday <- paste0(0,pdata$SamplingDay)
len <- nchar(sday)
sday <- substr(sday,len-1,len)
what <- paste(pdata$Treatment,sday,sep="")
X <- list(status= what)
names(X[[1]]) <- rownames(pdata)
X## $status
## C_S1_10 C_S1_7 C_S1_8 C_S1_9 C_S7_11 C_S7_12 C_S7_13
## "C01" "C01" "C01" "C01" "C07" "C07" "C07"
## C_S7_14 C_S8_15 C_S8_16 C_S8_17 C_S8_18 C_S14_19 C_S14_20
## "C07" "C08" "C08" "C08" "C08" "C14" "C14"
## C_S14_21 C_S14_22 H_S1_10 H_S1_7 H_S1_8 H_S1_9 H_S7_11
## "C14" "C14" "H01" "H01" "H01" "H01" "H07"
## H_S7_12 H_S7_13 H_S7_14 H_S8_15 H_S8_16 H_S8_17 H_S8_18
## "H07" "H07" "H07" "H08" "H08" "H08" "H08"
## H_S14_19 H_S14_20 H_S14_21 H_S14_22
## "H14" "H14" "H14" "H14"table(pdata$SamplingDay, what)## what
## C01 C07 C08 C14 H01 H07 H08 H14
## 1 4 0 0 0 4 0 0 0
## 7 0 4 0 0 0 4 0 0
## 8 0 0 4 0 0 0 4 0
## 14 0 0 0 4 0 0 0 4Put datasets into list X and ensure that they all have same order of samples, as in phenodata.
datanames## [1] "hormonomics" "metabolomics" "qPCR"for(i in 1:length(datanames)){
x <- get(datanames[i])
rownames(x) <- x[,1]
x <- x[,-1]
X[[i+1]] <- x[rownames(pdata),]
names(X)[i+1] <- datanames[i]
}
str(X)## List of 4
## $ status : Named chr [1:32] "C01" "C01" "C01" "C01" ...
## ..- attr(*, "names")= chr [1:32] "C_S1_10" "C_S1_7" "C_S1_8" "C_S1_9" ...
## $ hormonomics :'data.frame': 32 obs. of 12 variables:
## ..$ IAA : num [1:32] 37.2 45.9 47.6 37.6 49.5 ...
## ..$ oxIAA : num [1:32] 61.5 67.8 52.5 48.7 76.8 ...
## ..$ IAA.Asp : num [1:32] 2.22 1.99 1.5 2.24 1.93 ...
## ..$ ABA : num [1:32] 36.8 33.1 41.7 39.1 80 ...
## ..$ PA : num [1:32] 92 92.6 93.8 91.4 231.7 ...
## ..$ DPA : num [1:32] 45.6 62.1 55.1 59.8 178.1 ...
## ..$ SA : num [1:32] 505 519 275 628 1315 ...
## ..$ JA : num [1:32] 2.69 2.8 4.76 4.62 6.1 ...
## ..$ JA.Ile : num [1:32] 0.553 0.566 0.427 0.203 0.623 ...
## ..$ X9.10.dhJA: num [1:32] 5.45 3.7 2.88 5.17 5.01 ...
## ..$ X12.OH.JA : num [1:32] 16.1 24 27.9 25.7 244.1 ...
## ..$ cisOPDA : num [1:32] 354 403 645 750 1295 ...
## $ metabolomics:'data.frame': 32 obs. of 22 variables:
## ..$ Glukose : num [1:32] 2.13 2.2 0.82 2.55 4.77 6.3 7.24 3.09 6.34 9.49 ...
## ..$ Fructose: num [1:32] 2.7 2.9 1.59 3.01 3.97 7.16 5.41 4.04 8.52 7.75 ...
## ..$ Sucrose : num [1:32] 3.45 3.38 2.45 4.06 3.53 ...
## ..$ Starch : num [1:32] 22.06 12.74 9.98 15.13 16.25 ...
## ..$ Asp : num [1:32] 1040 844 887 988 793 ...
## ..$ Glu : num [1:32] 2514 1966 2068 2348 2109 ...
## ..$ Asn : num [1:32] 178 168 167 172 277 ...
## ..$ Ser : num [1:32] 598 498 441 538 368 ...
## ..$ Gln : num [1:32] 498 409 400 466 266 ...
## ..$ Gly : num [1:32] 137.8 104.6 92.1 117.2 68.7 ...
## ..$ His : num [1:32] 20.7 13.3 17.1 16.8 13.7 19.4 20.5 17.8 8.8 18.4 ...
## ..$ Arg : num [1:32] 26.9 29.1 29.5 24.9 38.3 63.9 38.6 47.3 54.7 67.2 ...
## ..$ Thr : num [1:32] 225 208 216 252 189 ...
## ..$ Ala : num [1:32] 702 496 515 653 296 ...
## ..$ Pro : num [1:32] 48.6 57.3 53.5 58.6 85.8 ...
## ..$ Tyr : num [1:32] 25.3 22.1 24.4 22.3 31.2 40.6 50.9 37.1 20.2 34.6 ...
## ..$ Val : num [1:32] 54.1 51.3 52.9 58.4 73.1 ...
## ..$ Met : num [1:32] 7.3 6.9 8.6 7.9 4.7 8.6 4.5 6.6 3.8 2.2 ...
## ..$ Ile : num [1:32] 48.8 47 54.3 52.9 60.4 80.6 62.7 63 29.9 48.3 ...
## ..$ Lys : num [1:32] 26 21.8 26.5 26.7 28.4 33.1 29.2 38.4 20.6 41.8 ...
## ..$ Leu : num [1:32] 18.8 15.7 19.5 20.1 27.4 27.1 21.3 26.3 38.5 45.7 ...
## ..$ Phe : num [1:32] 119 88 86.8 111 98.3 ...
## $ qPCR :'data.frame': 32 obs. of 14 variables:
## ..$ RbohA : num [1:32] 1.35 1.5 1.26 1.43 1.4 ...
## ..$ SnRK2 : num [1:32] 1.5 1.63 1.63 1.53 1.12 ...
## ..$ ACO2 : num [1:32] 0.15 0.196 0.44 0.177 0.537 ...
## ..$ HSP70 : num [1:32] 1.013 1.067 0.883 1.042 1.028 ...
## ..$ PR1b : num [1:32] 0.324 0.154 0.269 0.256 0.227 ...
## ..$ RD29B : num [1:32] 0.017 0.0465 0.017 0.0824 1.7519 ...
## ..$ X13.LOX: num [1:32] 0.7 0.66 0.766 0.734 1.136 ...
## ..$ P5CS : num [1:32] 3.648 3.265 2.152 3.156 0.863 ...
## ..$ ERF1 : num [1:32] 0.56 0.638 0.586 0.664 1.635 ...
## ..$ CAT1 : num [1:32] 0.64 0.635 0.687 0.705 0.811 ...
## ..$ CO : num [1:32] 2.42 5.56 2.74 2.47 1.71 ...
## ..$ SWEET : num [1:32] 0.816 1.874 0.931 0.934 1.495 ...
## ..$ SP6A : num [1:32] 0.122 0.239 0.33 0.122 3.386 ...
## ..$ M0ZJG3 : num [1:32] 1.211 1.376 1.007 0.903 2.358 ...names(X)## [1] "status" "hormonomics" "metabolomics" "qPCR"Same sample names?
OK <- TRUE
for(i in 2:length(X)) {
print(ok <- all(names(X[[1]])==rownames(X[[i]])))
OK <- OK&ok
}## [1] TRUE
## [1] TRUE
## [1] TRUESample names in datasets match.
Put data into safe object DATAall.
DATAall <- X
DATA <- XWe will also need the names of treatment groups.
groups <- unique(pdata$Treatment)
groups## [1] "C" "H"5 Multiomics data integration with DIABLO
CCDATA <- DATA
names(CCDATA)## [1] "status" "hormonomics" "metabolomics" "qPCR"write("Entering 035-DIABLO !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!", "bla.log", append=!TRUE)
write(.oroot, "bla.log", append=TRUE)
write(pisa$oroot, "bla.log", append=TRUE)
write(commandArgs(trailingOnly = TRUE), "bla.log",append=TRUE)out <- ""
out <- paste(out,knit_child("035-DIABLO-2.Rmd", quiet=TRUE))
cat(out)Child: 035-DIABLO-2.Rmd ## DIABLO hormonomics, metabolomics, qPCR ## DIABLO
DIABLO from mixOmics enables integration of more than two datasets.
5.1 Organize data
Thre datasets are organized as a list of matrices with same samples as rows and variables in columns.
data <- CCDATA[-1]
str(data)## List of 3
## $ hormonomics :'data.frame': 32 obs. of 12 variables:
## ..$ IAA : num [1:32] 37.2 45.9 47.6 37.6 49.5 ...
## ..$ oxIAA : num [1:32] 61.5 67.8 52.5 48.7 76.8 ...
## ..$ IAA.Asp : num [1:32] 2.22 1.99 1.5 2.24 1.93 ...
## ..$ ABA : num [1:32] 36.8 33.1 41.7 39.1 80 ...
## ..$ PA : num [1:32] 92 92.6 93.8 91.4 231.7 ...
## ..$ DPA : num [1:32] 45.6 62.1 55.1 59.8 178.1 ...
## ..$ SA : num [1:32] 505 519 275 628 1315 ...
## ..$ JA : num [1:32] 2.69 2.8 4.76 4.62 6.1 ...
## ..$ JA.Ile : num [1:32] 0.553 0.566 0.427 0.203 0.623 ...
## ..$ X9.10.dhJA: num [1:32] 5.45 3.7 2.88 5.17 5.01 ...
## ..$ X12.OH.JA : num [1:32] 16.1 24 27.9 25.7 244.1 ...
## ..$ cisOPDA : num [1:32] 354 403 645 750 1295 ...
## $ metabolomics:'data.frame': 32 obs. of 22 variables:
## ..$ Glukose : num [1:32] 2.13 2.2 0.82 2.55 4.77 6.3 7.24 3.09 6.34 9.49 ...
## ..$ Fructose: num [1:32] 2.7 2.9 1.59 3.01 3.97 7.16 5.41 4.04 8.52 7.75 ...
## ..$ Sucrose : num [1:32] 3.45 3.38 2.45 4.06 3.53 ...
## ..$ Starch : num [1:32] 22.06 12.74 9.98 15.13 16.25 ...
## ..$ Asp : num [1:32] 1040 844 887 988 793 ...
## ..$ Glu : num [1:32] 2514 1966 2068 2348 2109 ...
## ..$ Asn : num [1:32] 178 168 167 172 277 ...
## ..$ Ser : num [1:32] 598 498 441 538 368 ...
## ..$ Gln : num [1:32] 498 409 400 466 266 ...
## ..$ Gly : num [1:32] 137.8 104.6 92.1 117.2 68.7 ...
## ..$ His : num [1:32] 20.7 13.3 17.1 16.8 13.7 19.4 20.5 17.8 8.8 18.4 ...
## ..$ Arg : num [1:32] 26.9 29.1 29.5 24.9 38.3 63.9 38.6 47.3 54.7 67.2 ...
## ..$ Thr : num [1:32] 225 208 216 252 189 ...
## ..$ Ala : num [1:32] 702 496 515 653 296 ...
## ..$ Pro : num [1:32] 48.6 57.3 53.5 58.6 85.8 ...
## ..$ Tyr : num [1:32] 25.3 22.1 24.4 22.3 31.2 40.6 50.9 37.1 20.2 34.6 ...
## ..$ Val : num [1:32] 54.1 51.3 52.9 58.4 73.1 ...
## ..$ Met : num [1:32] 7.3 6.9 8.6 7.9 4.7 8.6 4.5 6.6 3.8 2.2 ...
## ..$ Ile : num [1:32] 48.8 47 54.3 52.9 60.4 80.6 62.7 63 29.9 48.3 ...
## ..$ Lys : num [1:32] 26 21.8 26.5 26.7 28.4 33.1 29.2 38.4 20.6 41.8 ...
## ..$ Leu : num [1:32] 18.8 15.7 19.5 20.1 27.4 27.1 21.3 26.3 38.5 45.7 ...
## ..$ Phe : num [1:32] 119 88 86.8 111 98.3 ...
## $ qPCR :'data.frame': 32 obs. of 14 variables:
## ..$ RbohA : num [1:32] 1.35 1.5 1.26 1.43 1.4 ...
## ..$ SnRK2 : num [1:32] 1.5 1.63 1.63 1.53 1.12 ...
## ..$ ACO2 : num [1:32] 0.15 0.196 0.44 0.177 0.537 ...
## ..$ HSP70 : num [1:32] 1.013 1.067 0.883 1.042 1.028 ...
## ..$ PR1b : num [1:32] 0.324 0.154 0.269 0.256 0.227 ...
## ..$ RD29B : num [1:32] 0.017 0.0465 0.017 0.0824 1.7519 ...
## ..$ X13.LOX: num [1:32] 0.7 0.66 0.766 0.734 1.136 ...
## ..$ P5CS : num [1:32] 3.648 3.265 2.152 3.156 0.863 ...
## ..$ ERF1 : num [1:32] 0.56 0.638 0.586 0.664 1.635 ...
## ..$ CAT1 : num [1:32] 0.64 0.635 0.687 0.705 0.811 ...
## ..$ CO : num [1:32] 2.42 5.56 2.74 2.47 1.71 ...
## ..$ SWEET : num [1:32] 0.816 1.874 0.931 0.934 1.495 ...
## ..$ SP6A : num [1:32] 0.122 0.239 0.33 0.122 3.386 ...
## ..$ M0ZJG3 : num [1:32] 1.211 1.376 1.007 0.903 2.358 ...length(data)## [1] 3In addition, outcome, phenotypic state or in our case treatment can also be determined. Let us say, the we have four treatments A, B, C and D.
state <- factor(CCDATA[[1]])
table(state)## state
## C01 C07 C08 C14 H01 H07 H08 H14
## 4 4 4 4 4 4 4 4str(state)## Factor w/ 8 levels "C01","C07","C08",..: 1 1 1 1 2 2 2 2 3 3 ...
## - attr(*, "names")= chr [1:32] "C_S1_10" "C_S1_7" "C_S1_8" "C_S1_9" ...5.2 Initial analysis
Note: Here we will follow the analysis as in http://mixomics.org/mixdiablo/diablo-tcga-case-study/.
5.3 Pairwise PLS comparisons
list.keepX = c(25, 25) # select arbitrary values of features to keep
list.keepY = c(25, 25)par(mfrow=c(2,2))
cutoffs <- c(0.5)
pairs <- combn(1:length(names(data)),2)
nms <- names(data)
outn <- ""
j <- 1
ncomp <- length(data)
cols <- c('orange1', 'lightgreen', "red")
if(length(data)==3) pick <- 1:3 else pick <- c(4,1:3)
cols <- c('orange1', 'brown1', 'lightgreen',"lightblue")[pick]
pchs <- c(16, 17, 15, 18)[pick]
j <- 4
cutoff <- 0.5
for(j in 1:ncol(pairs) ){
pair <- pairs[,j]
pair
X <- CCDATA[[pair[1]+1]]
Y <- CCDATA[[pair[2]+1]]
list.keepX <- rep(min(ncol(X), 25), ncomp)
list.keepY <- rep(min(ncol(Y), 25), ncomp)
x <- spls(X, Y, ncomp=ncomp, keepX = list.keepX, keepY = list.keepY)
assign(paste0("spls",j),x)
cat("\n",paste(nms[pair]), "\n")
cat("Results in:",paste0("spls",j),"\n")
cat("Correlation between pls variates:\n")
print(round(cor(x$variates$X, x$variates$Y),5))
# for(cutoff in cutoffs){
plotVar(x, cutoff = cutoff, title = paste(nms[pair],collapse=", "),
legend = c(nms[pair][1], nms[pair][2]),
var.names = FALSE, style = 'graphics',
pch = pchs[pair], cex = c(2,2),
col = cols[pair])
# }
}##
## hormonomics metabolomics
## Results in: spls1
## Correlation between pls variates:
## comp1 comp2 comp3
## comp1 0.91126 0.00000 0.00000
## comp2 0.13029 0.80764 0.00000
## comp3 0.00165 -0.00098 0.61826##
## hormonomics qPCR
## Results in: spls2
## Correlation between pls variates:
## comp1 comp2 comp3
## comp1 0.69330 0.00000 0.000
## comp2 -0.40812 0.63635 0.000
## comp3 -0.31053 0.39759 0.549##
## metabolomics qPCR
## Results in: spls3
## Correlation between pls variates:
## comp1 comp2 comp3
## comp1 0.74472 0.00000 0.00000
## comp2 0.38516 0.60146 0.00000
## comp3 0.27693 0.40129 0.61882
Circle Correlation Plots for pairwise PLS models on ADAPT data. Only displays the at most top 25 features for each dimension, subsetting by those with a correlation above 0.5.
5.4 Initial DIABLO model
Following the suggestion in the source, we will use design matrices with values. This is supposed to keep low classification error rate.
entry <- .entry
design = matrix(entry, ncol = length(data), nrow = length(data),
dimnames = list(names(data), names(data)))
diag(design) = 0 # set diagonal to 0s
design## hormonomics metabolomics qPCR
## hormonomics 0.0 0.5 0.5
## metabolomics 0.5 0.0 0.5
## qPCR 0.5 0.5 0.0With a design in place, the initial DIABLO model can be generated. An arbitrarily high number of components (ncomp = 5) will be used.
# form basic DIABLO model
Y <- state
basic.diablo.model = block.splsda(X = data, Y = Y, ncomp = 5, design = design)## Design matrix has changed to include Y; each block will be
## linked to Y.5.5 Tuning the number of components
To choose the number of components for the final DIABLO model, the function perf() is run with 3-fold cross-validation repeated 10 times. Fold number should be smaller than minimal number of samples in groups.
# run component number tuning with repeated CV
system.time(perf.diablo = perf(basic.diablo.model, validation = 'Mfold',
folds = 3, nrepeat = 10))
plot(perf.diablo) # plot output of tuningThe classification rate stays
# set the optimal ncomp value
ncomp <- perf.diablo$choice.ncomp$WeightedVote["Overall.BER", "centroids.dist"]
# show the optimal choice for ncomp for each dist metric
perf.diablo$choice.ncomp$WeightedVoteFor classification, the analysis suggests the number of components.
(Tuning takes some time.) From previous tuning sessions one can conclude, that the classification rate stays roughly unchanged after two to four components, so we will set the number of components to number of data sets:
ncomp <- length(data)
ncomp## [1] 35.6 Tuning the number of features
This tuning function should be used to tune the keepX parameters in the block.splsda() function.
We choose the optimal number of variables to select in each data set using the tune.block.splsda() function, for a grid of keepX values for each type of omics. Note that the function has been set to favour a relatively small signature while allowing us to obtain a sufficient number of variables for downstream validation and/or interpretation. See ?tune.block.splsda.
The function tune is run with 10-fold cross validation, but repeated only once. Note that for a more thorough tuning process, provided sufficient computational time, we could increase the nrepeat argument. Here we have saved the results into an RData object that is available for download as the tuning can take a very long time, especially on lower end machines.
x <- list()
for (i in 1:length(data)){
x[[i]] <- c( seq(5,min(30, ncol(data[[i]])) ,5))
}
names(x) <- names(data)
test.keepX <- x
test.keepX## $hormonomics
## [1] 5 10
##
## $metabolomics
## [1] 5 10 15 20
##
## $qPCR
## [1] 5 10#list (c(5:9, seq(10, 18, 2), seq(20,30,5)),
# c(5:9, seq(10, 18, 2), seq(20,30,5)),
# c(5:9, seq(10, 18, 2), seq(20,30,5)))# run the feature selection tuning
system.time(tune.model <- tune.block.splsda(X = data, Y = Y, ncomp = ncomp, cpus=4,
test.keepX = test.keepX, design = design,
validation = 'Mfold', folds = 3, nrepeat = 1,
dist = "centroids.dist")
)# run the feature selection tuning
system.time(tune.model <- tune.block.splsda(X = data, Y = Y, ncomp = ncomp, cpus=4,
test.keepX = test.keepX, design = design,
validation = 'loo', folds = 3, nrepeat = 1,
dist = "centroids.dist")
)The number of features to select on each component is returned in
list.keepX = tune.model$choice.keepX # set the optimal values of features to retain
list.keepXTuning of features can take a substantial amount of time. From previous runs, we can conclude to limit number of variates mostly in proteomics. However, I do not see the effect of this limitations except in a short list of, presumably, important features. Loadings are different. Previous analyses suggest the following list:
$metabolomics [1] 10 10 5 5 $hormonomics [1] 5 5 10 5 $qPCR [1] 10 10 5 5
keepX <- list(
metabolomics = rep(10, ncomp),
hormonomics = rep(10, ncomp),
qPCR = rep(10, ncomp)
)
list.keepX = list()
for (i in 1:length(data)) list.keepX[[i]] <- keepX[[names(data)[i]]]
names(list.keepX) <- names(data)
list.keepX## $hormonomics
## [1] 10 10 10
##
## $metabolomics
## [1] 10 10 10
##
## $qPCR
## [1] 10 10 105.7 Final model
The final DIABLO model is run as:
# set the optimised DIABLO model
final.diablo.model = block.splsda(X = data, Y = Y, ncomp = ncomp
, keepX = list.keepX
, design = design)## Design matrix has changed to include Y; each block will be
## linked to Y.The selected variables can be extracted with the function selectVar(), for example in each block, as seen below. Note that the stability of selected variables can be extracted from the output of the perf() function.
# the features selected from components
for (comp in 1:ncomp){
cat("\nComponent ", comp,":\n")
for(i in 1:length(data)){
cat(names(data)[i],"\n")
print(selectVar(final.diablo.model, comp = comp)[[i]]$name)
}
}##
## Component 1 :
## hormonomics
## [1] "DPA" "SA" "PA" "X12.OH.JA" "X9.10.dhJA"
## [6] "JA.Ile" "ABA" "IAA" "IAA.Asp" "cisOPDA"
## metabolomics
## [1] "Glukose" "Fructose" "Val" "Ile" "Tyr" "Lys"
## [7] "His" "Gln" "Pro" "Met"
## qPCR
## [1] "X13.LOX" "PR1b" "CAT1" "SP6A" "M0ZJG3" "HSP70"
## [7] "SWEET" "RbohA" "SnRK2" "ERF1"
##
## Component 2 :
## hormonomics
## [1] "ABA" "oxIAA" "IAA" "cisOPDA" "JA.Ile"
## [6] "PA" "IAA.Asp" "SA" "DPA" "X12.OH.JA"
## metabolomics
## [1] "Starch" "Ser" "Asn" "His" "Met" "Arg"
## [7] "Sucrose" "Gly" "Pro" "Glukose"
## qPCR
## [1] "SP6A" "SnRK2" "RD29B" "CO" "P5CS" "PR1b" "HSP70"
## [8] "M0ZJG3" "ACO2" "RbohA"
##
## Component 3 :
## hormonomics
## [1] "JA.Ile" "JA" "oxIAA" "IAA.Asp" "X9.10.dhJA"
## [6] "DPA" "PA" "SA" "IAA" "cisOPDA"
## metabolomics
## [1] "Met" "Gln" "Ala" "Leu" "Phe" "Gly" "Glu" "Ile" "Arg" "Asp"
## qPCR
## [1] "ACO2" "PR1b" "X13.LOX" "RD29B" "CO" "M0ZJG3"
## [7] "SnRK2" "P5CS" "ERF1" "SWEET"5.8 Sample plots
plotDIABLO() is a diagnostic plot to check whether the correlation between components from each data set has been maximised as specified in the design matrix. We specify which dimension to be assessed with the ncomp argument.
for(comp in 1:ncomp){
plotDiablo(final.diablo.model, ncomp = comp)
title(paste("Component",comp), adj=0.1, line=-1, outer=TRUE)
}
The sample plot with the plotIndiv() function projects each sample into the space spanned by the components of each block. Clustering of the samples can be assessed with this plot.
plind <- plotIndiv(final.diablo.model, ind.names = FALSE, legend = TRUE,
title = 'DIABLO Sample Plots'
, ellipse = TRUE
)## Warning: It is deprecated to specify `guide = FALSE` to remove a
## guide. Please use `guide = "none"` instead.
In the arrow plot below (Figure), the start of the arrow indicates the centroid between all data sets for a given sample and the tips of the arrows indicate the location of that sample in each block. Such graphics highlight the agreement between all data sets at the sample level. While somewhat difficult to interpret, even qualitatively, Figure 5 shows that the agreement within the LumA group seems to be the highest and lowest in the Her2 group.
plotArrow(final.diablo.model, ind.names = FALSE, legend = TRUE,
title = paste(groups,collapse=", ")
)
5.9 Variable plots
Several graphical outputs are available to visualise and mine the associations between the selected variables.
The best starting point to evaluate the correlation structure between variables is with the correlation circle plot, depicted in Figure 6. A majority of the miRNA variables are positively correlated with the first component while the mRNA variables seem to separate along this dimension. These first two components correlate highly with the selected variables from the proteomics dataset. From this, the correlation of each selected feature from all three datasets can be evaluated based on their proximity. s
if(length(data)==3) pick <- 1:3 else pick <- c(4,1:3)
cols <- c('orange1', 'brown1', 'lightgreen',"lightblue")[pick]
pchs <- c(16, 17, 15, 18)[pick]
plotVar(final.diablo.model, var.names = FALSE,
style = 'graphics', legend = TRUE
, pch = pchs, cex = rep(2,length(data))
, col = cols
)
The circos plot is exclusive to integrative frameworks and represents the correlations between variables of different types, represented on the side quadrants. From Figure 7, it seems that the miRNA variables are almost entirely negatively correlated with the other two dataframes. The proteomics features are the opposite, such that they display primarily positive correlations while the mRNA variables are more mixed. Note that these correlations are above a value of 0.7 (cutoff = 0.7). All the interpretations made above are only relevant for features with very strong correlations.
circosPlot(final.diablo.model, cutoff = 0.7, line = TRUE,
color.blocks= cols,
color.cor = c(3,2), size.labels = 1
, xpd=TRUE)
5.10 Relevance networks
Another visualisation of the correlations between the different types of variables is the relevance network, which is also built on the similarity matrix (as is the circos plot). Each colour represents a type of variable.
blocks <- combn(length(data),2)
j <- 1
cutoff <- 0.8
out35a <- ""
for(j in 1:ncol(blocks)){
out35a <- paste( out35a, knit_child("035a-DIABLO-network.Rmd", quiet=!TRUE))
if(interactive()) readline()
}
cat(out35a)5.10.1 hormonomics and metabolomics
nfn <- paste0("network-035a-",paste(names(data)[blocks[,j]], collapse="-"),"-",cutoff*10)
#nfn <- paste0("network-035a-",j,"-",cutoff*10)
write(nfn, "bla.log", append=TRUE)Error: figure margins too large
png(paste0(nfn,".png"), res = 600, width = 4000, height = 4000)
nw <- network(final.diablo.model
, blocks = blocks[,j]
, color.node = cols[blocks[,j]]
, cutoff = cutoff
, shape.node = "rectangle"
, save = "png"
, name.save = nfn
)
#title(main=paste(names(data)[blocks[,j]], sep=", "),
#sub=paste("Cutoff = ",cutoff))
#
#dev.off()Image: 
5.10.2 hormonomics and qPCR
nfn <- paste0("network-035a-",paste(names(data)[blocks[,j]], collapse="-"),"-",cutoff*10)
#nfn <- paste0("network-035a-",j,"-",cutoff*10)
write(nfn, "bla.log", append=TRUE)Error: figure margins too large
png(paste0(nfn,".png"), res = 600, width = 4000, height = 4000)
nw <- network(final.diablo.model
, blocks = blocks[,j]
, color.node = cols[blocks[,j]]
, cutoff = cutoff
, shape.node = "rectangle"
, save = "png"
, name.save = nfn
)
#title(main=paste(names(data)[blocks[,j]], sep=", "),
#sub=paste("Cutoff = ",cutoff))
#
#dev.off()Image: 
5.10.3 metabolomics and qPCR
nfn <- paste0("network-035a-",paste(names(data)[blocks[,j]], collapse="-"),"-",cutoff*10)
#nfn <- paste0("network-035a-",j,"-",cutoff*10)
write(nfn, "bla.log", append=TRUE)Error: figure margins too large
png(paste0(nfn,".png"), res = 600, width = 4000, height = 4000)
nw <- network(final.diablo.model
, blocks = blocks[,j]
, color.node = cols[blocks[,j]]
, cutoff = cutoff
, shape.node = "rectangle"
, save = "png"
, name.save = nfn
)
#title(main=paste(names(data)[blocks[,j]], sep=", "),
#sub=paste("Cutoff = ",cutoff))
#
#dev.off()Image: 
The function plotLoadings() visualises the loading weights of each selected variable on each component and each data set (Figure 9). The colour indicates the class in which the variable has the maximum level of expression (contrib = ‘max’) using the median (method = ‘median’). Figure 9 depicts the loading values for the second dimension.
5.11 Multipartite network
cutoff <- 0.0
x <- final.diablo.model
layout.fun <- NULL
label <- paste(.treat, collapse=", ")
out35b <- ""
out35b <- paste( out35b, knit_child("035b-multipartite-network.Rmd", quiet=TRUE))
cat(out35b)5.11.1 Cutoff = 0
ndata <- length(data)
lbl <- gsub(", ","-",label)
nfn <- paste("network-035b",lbl,cutoff*10,sep="-")
#png(nfn, res = 600, width = 4000, height = 4000)
write(nfn, "bla.log", append=TRUE)set.seed(1234)
nw <- my.network(x
, blocks = 1:ndata
, color.node = cols
, cutoff = cutoff
, shape.node = "rectangle"
, layout = layout.fun
, save = "png"
, name.save = nfn
)
# title( #main=paste(names(data), sep=", "),
# sub=paste("Cutoff = ",cutoff))
# title(label,adj=0.8,outer=TRUE,line=-1)
# legend("bottomright", pch=15,pt.cex=2,col=cols, legend=names(data),
# bty="n")
# text(ly[,1],ly[,2],names(V(nw$gR)))
#dev.off()
network-035b-C-H-0
if(exists(deparse(substitute(nw)))) ly <- nw$layout else ly <- NULLcutoff <- 0.5
x <- final.diablo.model
layout.fun <- my.layout
label <- paste(.treat, collapse=", ")
out35b <- ""
out35b <- paste( out35b, knit_child("035b-multipartite-network.Rmd", quiet=TRUE))
cat(out35b)5.11.2 Cutoff = 0.5
ndata <- length(data)
lbl <- gsub(", ","-",label)
nfn <- paste("network-035b",lbl,cutoff*10,sep="-")
#png(nfn, res = 600, width = 4000, height = 4000)
write(nfn, "bla.log", append=TRUE)set.seed(1234)
nw <- my.network(x
, blocks = 1:ndata
, color.node = cols
, cutoff = cutoff
, shape.node = "rectangle"
, layout = layout.fun
, save = "png"
, name.save = nfn
)
# title( #main=paste(names(data), sep=", "),
# sub=paste("Cutoff = ",cutoff))
# title(label,adj=0.8,outer=TRUE,line=-1)
# legend("bottomright", pch=15,pt.cex=2,col=cols, legend=names(data),
# bty="n")
# text(ly[,1],ly[,2],names(V(nw$gR)))
#dev.off()
network-035b-C-H-5
5.11.3 Cutoff = 0.5
ndata <- length(data)
lbl <- gsub(", ","-",label)
nfn <- paste("network-035b",lbl,cutoff*10,sep="-")
#png(nfn, res = 600, width = 4000, height = 4000)
write(nfn, "bla.log", append=TRUE)set.seed(1234)
nw <- my.network(x
, blocks = 1:ndata
, color.node = cols
, cutoff = cutoff
, shape.node = "rectangle"
, layout = layout.fun
, save = "png"
, name.save = nfn
)
# title( #main=paste(names(data), sep=", "),
# sub=paste("Cutoff = ",cutoff))
# title(label,adj=0.8,outer=TRUE,line=-1)
# legend("bottomright", pch=15,pt.cex=2,col=cols, legend=names(data),
# bty="n")
# text(ly[,1],ly[,2],names(V(nw$gR)))
#dev.off()network-035b-C-H-5
cutoff <- 0.7
x <- final.diablo.model
label <- paste(.treat, collapse=", ")
out35b <- ""
out35b <- paste( out35b, knit_child("035b-multipartite-network.Rmd", quiet=TRUE))
cat(out35b)5.11.4 Cutoff = 0.7
ndata <- length(data)
lbl <- gsub(", ","-",label)
nfn <- paste("network-035b",lbl,cutoff*10,sep="-")
#png(nfn, res = 600, width = 4000, height = 4000)
write(nfn, "bla.log", append=TRUE)set.seed(1234)
nw <- my.network(x
, blocks = 1:ndata
, color.node = cols
, cutoff = cutoff
, shape.node = "rectangle"
, layout = layout.fun
, save = "png"
, name.save = nfn
)
# title( #main=paste(names(data), sep=", "),
# sub=paste("Cutoff = ",cutoff))
# title(label,adj=0.8,outer=TRUE,line=-1)
# legend("bottomright", pch=15,pt.cex=2,col=cols, legend=names(data),
# bty="n")
# text(ly[,1],ly[,2],names(V(nw$gR)))
#dev.off()
network-035b-C-H-7
cutoff <- 0.8
x <- final.diablo.model
label <- paste(.treat, collapse=", ")
out35b <- ""
out35b <- paste( out35b, knit_child("035b-multipartite-network.Rmd", quiet=TRUE))
cat(out35b)5.11.5 Cutoff = 0.8
ndata <- length(data)
lbl <- gsub(", ","-",label)
nfn <- paste("network-035b",lbl,cutoff*10,sep="-")
#png(nfn, res = 600, width = 4000, height = 4000)
write(nfn, "bla.log", append=TRUE)set.seed(1234)
nw <- my.network(x
, blocks = 1:ndata
, color.node = cols
, cutoff = cutoff
, shape.node = "rectangle"
, layout = layout.fun
, save = "png"
, name.save = nfn
)
# title( #main=paste(names(data), sep=", "),
# sub=paste("Cutoff = ",cutoff))
# title(label,adj=0.8,outer=TRUE,line=-1)
# legend("bottomright", pch=15,pt.cex=2,col=cols, legend=names(data),
# bty="n")
# text(ly[,1],ly[,2],names(V(nw$gR)))
#dev.off()
network-035b-C-H-8
for(i in 1:ncomp)
plotLoadings(final.diablo.model, comp = i, contrib = 'max', method = 'median')
## Heatmap (clustered image map)
The cimDIABLO() function is a clustered image map specifically implemented to represent the multi-omics molecular signature expression for each sample. From figure below the areas of homogeneous expression levels for a set of samples across a set of features can be determined. For instance, the H14 samples were the only group to show extremely high levels of expression for a specific set of genes and metabolites. This indicates these features are fairly discriminating for this subtype.
cimfn <- "cim.png"
png(cimfn, res = 600, width = 4000, height = 4000)
cimDiablo(final.diablo.model, size.legend=0.7)
dev.off()## pdf
## 2
cim.png
An AUC plot per block can also be obtained using the function auroc(). The interpretation of this output may not be particularly insightful in relation to the performance evaluation of our methods, but can complement the statistical analysis..
par(mfrow=c(2,2))
for(i in 1:length(data))
auc.splsda = auroc(final.diablo.model, roc.block = names(data[i]),
roc.comp = 2, print = FALSE)
res12 <- final.diablo.model5.12 Final analysis
Design matrix determines which datasets (blocks) are connected. We will connect all blocks:
According to recommendation from tuning, design matrix should have small antries that provide better predictions.
entry <- .entry
design <- matrix(entry, length(data), length(data),
dimnames = list(names(data), names(data)))
diag(design) <- 0
design## hormonomics metabolomics qPCR
## hormonomics 0.0 0.5 0.5
## metabolomics 0.5 0.0 0.5
## qPCR 0.5 0.5 0.0Based on tuning part, we can keep two to four components.
list.keepX## $hormonomics
## [1] 10 10 10
##
## $metabolomics
## [1] 10 10 10
##
## $qPCR
## [1] 10 10 10ncomp = length(list.keepX[[1]])
ncomp## [1] 3You can also control the number of variables to be kept for each component.
#keepX <- list(XX=rep(3,ncomp), YY=rep(3,ncomp), ZZ=rep(3, ncomp))
#keepXWe will not restrict variables.
Now we can do the calculation
res <- block.splsda(X = data
, Y = as.factor(state)
, ncomp = ncomp
, keepX = list.keepX
, design = design
)## Design matrix has changed to include Y; each block will be
## linked to Y.Model res and final.diablo.model do not differ:
sapply(names(res$loadings),
function(name) all(res$loadings[[name]] -
final.diablo.model$loadings[[name]]==0))## hormonomics metabolomics qPCR Y
## TRUE TRUE TRUE TRUEEstimate classification error rate. The error rate should drop by more components used.
# run component number tuning with repeated CV
system.time(perf.diablo <- perf(res, validation = 'Mfold',
folds = 3, nrepeat = 10))
plot(perf.diablo) # plot output of tuningNames of kept variables
# the features selected to form components
for (comp in 1:ncomp){
cat("\nComponent ", comp,":\n")
for(i in 1:length(data)){
cat(names(data)[i],"\n")
print(selectVar(res, comp = comp)[[i]]$name)
}
}##
## Component 1 :
## hormonomics
## [1] "DPA" "SA" "PA" "X12.OH.JA" "X9.10.dhJA"
## [6] "JA.Ile" "ABA" "IAA" "IAA.Asp" "cisOPDA"
## metabolomics
## [1] "Glukose" "Fructose" "Val" "Ile" "Tyr" "Lys"
## [7] "His" "Gln" "Pro" "Met"
## qPCR
## [1] "X13.LOX" "PR1b" "CAT1" "SP6A" "M0ZJG3" "HSP70"
## [7] "SWEET" "RbohA" "SnRK2" "ERF1"
##
## Component 2 :
## hormonomics
## [1] "ABA" "oxIAA" "IAA" "cisOPDA" "JA.Ile"
## [6] "PA" "IAA.Asp" "SA" "DPA" "X12.OH.JA"
## metabolomics
## [1] "Starch" "Ser" "Asn" "His" "Met" "Arg"
## [7] "Sucrose" "Gly" "Pro" "Glukose"
## qPCR
## [1] "SP6A" "SnRK2" "RD29B" "CO" "P5CS" "PR1b" "HSP70"
## [8] "M0ZJG3" "ACO2" "RbohA"
##
## Component 3 :
## hormonomics
## [1] "JA.Ile" "JA" "oxIAA" "IAA.Asp" "X9.10.dhJA"
## [6] "DPA" "PA" "SA" "IAA" "cisOPDA"
## metabolomics
## [1] "Met" "Gln" "Ala" "Leu" "Phe" "Gly" "Glu" "Ile" "Arg" "Asp"
## qPCR
## [1] "ACO2" "PR1b" "X13.LOX" "RD29B" "CO" "M0ZJG3"
## [7] "SnRK2" "P5CS" "ERF1" "SWEET"One would like to reduce the number of nodes, especially for proteomics data. One option is to reduce datasets in a way to keep only the variables in the selectVars in original data in . We will keep variables from the first two components.
keptVars <- unique(c(
selectVar(res, comp=1)[[1]]$name
,selectVar(res, comp=2)[[1]]$name
)
)
which(keptVars%in%selectVar(res, comp=1)[[1]]$name)## [1] 1 2 3 4 5 6 7 8 9 10which(keptVars%in%selectVar(res, comp=2)[[1]]$name)## [1] 1 2 3 4 6 7 8 9 10 115.12.1 Visualise variables
Loadings
sapply(res$loadings, head, 30)## $hormonomics
## comp1 comp2 comp3
## IAA 0.036703482 -0.40678302 0.02662456
## oxIAA 0.000000000 0.41472146 -0.46097105
## IAA.Asp -0.007071215 -0.19836112 -0.38748591
## ABA 0.216332456 0.49851046 0.00000000
## PA 0.426177444 0.25246990 0.07019303
## DPA 0.549700457 -0.08274318 0.14734983
## SA 0.428211919 -0.17502093 0.03899104
## JA 0.000000000 0.00000000 0.48879193
## JA.Ile 0.246436078 -0.32453847 -0.55019009
## X9.10.dhJA 0.297909166 0.00000000 -0.25753363
## X12.OH.JA 0.367626768 0.05459029 0.00000000
## cisOPDA 0.003133045 0.40638498 -0.02327900
##
## $metabolomics
## comp1 comp2 comp3
## Glukose 0.48451843 0.0007744579 0.00000000
## Fructose 0.44913578 0.0000000000 0.00000000
## Sucrose 0.00000000 0.1755425740 0.00000000
## Starch 0.00000000 0.6230690186 0.00000000
## Asp 0.00000000 0.0000000000 0.02433324
## Glu 0.00000000 0.0000000000 -0.08664505
## Asn 0.00000000 -0.3852103155 0.00000000
## Ser 0.00000000 0.4187171275 0.00000000
## Gln -0.13859433 0.0000000000 0.52431197
## Gly 0.00000000 0.0875580979 0.12358545
## His 0.21679051 -0.3308628088 0.00000000
## Arg 0.00000000 -0.2436655212 -0.02652000
## Thr 0.00000000 0.0000000000 0.00000000
## Ala 0.00000000 0.0000000000 0.41741363
## Pro 0.11913529 0.0723114386 0.00000000
## Tyr 0.34922488 0.0000000000 0.00000000
## Val 0.37933829 0.0000000000 0.00000000
## Met -0.07239333 -0.2748085404 0.62202620
## Ile 0.35208370 0.0000000000 0.08509517
## Lys 0.29674894 0.0000000000 0.00000000
## Leu 0.00000000 0.0000000000 -0.27616600
## Phe 0.00000000 0.0000000000 0.23740522
##
## $qPCR
## comp1 comp2 comp3
## RbohA 0.16466558 -0.01682425 0.000000000
## SnRK2 0.13112994 -0.30963253 -0.139934444
## ACO2 0.00000000 -0.04461440 -0.883294093
## HSP70 0.24795487 -0.11801339 0.000000000
## PR1b 0.37477829 0.12308322 -0.306943422
## RD29B 0.00000000 0.28522220 0.155469415
## X13.LOX 0.58349667 0.00000000 0.168233520
## P5CS 0.00000000 -0.14974981 0.066857223
## ERF1 0.09174938 0.00000000 -0.065022964
## CAT1 0.37449858 0.00000000 0.000000000
## CO 0.00000000 -0.15196379 0.150498594
## SWEET 0.21148675 0.00000000 0.004392953
## SP6A 0.34570655 0.85834408 0.000000000
## M0ZJG3 0.31681955 -0.09567284 0.148846821
##
## $Y
## comp1 comp2 comp3
## C01 -0.346770705 0.03865545 0.25532621
## C07 -0.110583348 0.11426169 0.19174216
## C08 -0.008904022 0.38818375 0.02515787
## C14 0.063004541 0.64332045 -0.38125662
## H01 -0.287945254 -0.14936480 0.29508918
## H07 -0.116690109 -0.40997461 0.17110929
## H08 -0.065417492 -0.44511069 -0.76872677
## H14 0.873306389 -0.17997122 0.21155868#plotLoadings(res, comp = 1, method = 'median')
#plotLoadings(res, comp = 1, method = 'median', contrib="max")
for( i in 1:ncomp)
plotLoadings(res, comp = i, method = 'median', contrib="max")
Plot variables
#plotVar(res, cutoff=0.5, legend = TRUE, overlap=!FALSE, style='graphics')
#plotVar(res, cutoff=0.5, legend = TRUE, overlap=FALSE, style='graphics')
plotVar(res, cutoff=0.5, legend = TRUE, comp=c(1,2), overlap=FALSE, style='ggplot2', col=cols)
plotVar(res, cutoff=0.5, legend = TRUE, comp=c(2,3), overlap=FALSE, col=cols)
5.13 Differential networks
Here we will show differential networks between treatments.
cutoffs <- c(0.7)
pairs <- combn(1:length(names(res$X)),2)
outn <- ""
j <- 4
cutoff <- 0.5
for(j in 1:ncol(pairs) ){
pair <- pairs[,j]
X <- data[[pair[1]]]
Y <- data[[pair[2]]]
datasets <- names(data)[pair]
outn <- paste( outn, knit_child("023-prepare-networkdiff.Rmd", quiet=TRUE))
for(cutoff in cutoffs){
outn <- paste( outn, knit_child("035-Network.Rmd", quiet=TRUE))
}
}
cat(outn)size.variables <- 1
sim <-circosPlot(final.diablo.model, cutoff = 0.5, line = TRUE,
color.blocks= cols,
color.cor = c(3,2), size.labels = 1
, size.variables = size.variables
, xpd=TRUE)
circosPlot(final.diablo.model, cutoff = 0.78, line = TRUE,
color.blocks= cols,
color.cor = c(3,2), size.labels = 1
, size.variables = size.variables
, xpd=TRUE)
circosPlot(final.diablo.model, cutoff = 0.9, line = TRUE,
color.blocks= cols,
color.cor = c(3,2), size.labels = 1
, size.variables = size.variables
, xpd=TRUE)
circosPlot(final.diablo.model, cutoff = 0.95, line = TRUE,
color.blocks= cols,
color.cor = c(3,2), size.labels = 1
, size.variables = size.variables
, xpd=TRUE)
6 Export networks to file
Partial models for each treatment
6.1 Network for C.
filter <- pdata$Treatment %in% .treat[1]
XX1 <- lapply(CCDATA, function(x) if(is.null(dim(x))) x[filter] else x[filter,])
table(XX1$status)##
## C01 C07 C08 C14
## 4 4 4 4res1 <- block.splsda(X = XX1[-1]
, Y = as.factor(XX1[[1]])
, ncomp = ncomp
, keepX = list.keepX
, design = design
)## Design matrix has changed to include Y; each block will be
## linked to Y.cutoff <- 0.0
x <- res1
layout.fun <- NULL
label <-.treat[1]
out23b <- ""
out23b <- paste( out23b, knit_child("035b-multipartite-network.Rmd", quiet=TRUE))
N1 <- nw
cat(out23b)6.1.1 Cutoff = 0
ndata <- length(data)
lbl <- gsub(", ","-",label)
nfn <- paste("network-035b",lbl,cutoff*10,sep="-")
#png(nfn, res = 600, width = 4000, height = 4000)
write(nfn, "bla.log", append=TRUE)set.seed(1234)
nw <- my.network(x
, blocks = 1:ndata
, color.node = cols
, cutoff = cutoff
, shape.node = "rectangle"
, layout = layout.fun
, save = "png"
, name.save = nfn
)
# title( #main=paste(names(data), sep=", "),
# sub=paste("Cutoff = ",cutoff))
# title(label,adj=0.8,outer=TRUE,line=-1)
# legend("bottomright", pch=15,pt.cex=2,col=cols, legend=names(data),
# bty="n")
# text(ly[,1],ly[,2],names(V(nw$gR)))
#dev.off()
Save network layout for further plots, used by layout function my.layout.
ly <- nw$layoutcutoff <- 0.7
x <- res1
layout.fun <- my.layout
label <- .treat[1]
out23b <- ""
out23b <- paste( out23b, knit_child("035b-multipartite-network.Rmd", quiet=TRUE))
cat(out23b)6.1.2 Cutoff = 0.7
ndata <- length(data)
lbl <- gsub(", ","-",label)
nfn <- paste("network-035b",lbl,cutoff*10,sep="-")
#png(nfn, res = 600, width = 4000, height = 4000)
write(nfn, "bla.log", append=TRUE)set.seed(1234)
nw <- my.network(x
, blocks = 1:ndata
, color.node = cols
, cutoff = cutoff
, shape.node = "rectangle"
, layout = layout.fun
, save = "png"
, name.save = nfn
)
# title( #main=paste(names(data), sep=", "),
# sub=paste("Cutoff = ",cutoff))
# title(label,adj=0.8,outer=TRUE,line=-1)
# legend("bottomright", pch=15,pt.cex=2,col=cols, legend=names(data),
# bty="n")
# text(ly[,1],ly[,2],names(V(nw$gR)))
#dev.off()
network-035b-C-7
6.2 Network for H.
filter <- pdata$Treatment %in% .treat[2]
XX2 <- lapply(CCDATA, function(x) if(is.null(dim(x))) x[filter] else x[filter,])
table(XX2$status)##
## H01 H07 H08 H14
## 4 4 4 4res2 <- block.splsda(X = XX2[-1]
, Y = as.factor(XX2[[1]])
, ncomp = ncomp
, keepX = list.keepX
, design = design
)## Design matrix has changed to include Y; each block will be
## linked to Y.cutoff <- 0.0
x <- res2
layout.fun <- NULL
label <- .treat[2]
out23b <- ""
out23b <- paste( out23b, knit_child("035b-multipartite-network.Rmd", quiet=TRUE))
N2 <- nw
cat(out23b)6.2.1 Cutoff = 0
ndata <- length(data)
lbl <- gsub(", ","-",label)
nfn <- paste("network-035b",lbl,cutoff*10,sep="-")
#png(nfn, res = 600, width = 4000, height = 4000)
write(nfn, "bla.log", append=TRUE)set.seed(1234)
nw <- my.network(x
, blocks = 1:ndata
, color.node = cols
, cutoff = cutoff
, shape.node = "rectangle"
, layout = layout.fun
, save = "png"
, name.save = nfn
)
# title( #main=paste(names(data), sep=", "),
# sub=paste("Cutoff = ",cutoff))
# title(label,adj=0.8,outer=TRUE,line=-1)
# legend("bottomright", pch=15,pt.cex=2,col=cols, legend=names(data),
# bty="n")
# text(ly[,1],ly[,2],names(V(nw$gR)))
#dev.off()
network-035b-H-0
Save layout for further plots, used by layout function my.layout.
ly <- nw$layoutcutoff <- 0.7
x <- res2
layout.fun <- my.layout
label <- .treat[2]
out23b <- ""
out23b <- paste( out23b, knit_child("035b-multipartite-network.Rmd", quiet=TRUE))
cat(out23b)6.2.2 Cutoff = 0.7
ndata <- length(data)
lbl <- gsub(", ","-",label)
nfn <- paste("network-035b",lbl,cutoff*10,sep="-")
#png(nfn, res = 600, width = 4000, height = 4000)
write(nfn, "bla.log", append=TRUE)set.seed(1234)
nw <- my.network(x
, blocks = 1:ndata
, color.node = cols
, cutoff = cutoff
, shape.node = "rectangle"
, layout = layout.fun
, save = "png"
, name.save = nfn
)
# title( #main=paste(names(data), sep=", "),
# sub=paste("Cutoff = ",cutoff))
# title(label,adj=0.8,outer=TRUE,line=-1)
# legend("bottomright", pch=15,pt.cex=2,col=cols, legend=names(data),
# bty="n")
# text(ly[,1],ly[,2],names(V(nw$gR)))
#dev.off()
network-035b-H-7
Save network file for combined and single treatments. Networks are in objects res, res1 and res2.
6.3 Network for C and H.
# Complete network, cutoff = 0, both
datasets <- names(CCDATA[-1])
ndatasets<- length(datasets)
#
N12 <- network(res
, cutoff = 0
, blocks = 1:ndatasets
, shape.node = c("rectangle")
, save = "png"
, name.save="network-CH"
)
#
e <- extractEdges2(N12)
colnames(e)[ncol(e)] <- paste(.treat, collapse=".")
head(e)## edge group1 from
## ho.IAA_me.Glukose ho.IAA_me.Glukose hormonomics IAA
## ho.oxIAA_me.Glukose ho.oxIAA_me.Glukose hormonomics oxIAA
## ho.IAA.Asp_me.Glukose ho.IAA.Asp_me.Glukose hormonomics IAA.Asp
## ho.ABA_me.Glukose ho.ABA_me.Glukose hormonomics ABA
## ho.PA_me.Glukose ho.PA_me.Glukose hormonomics PA
## ho.DPA_me.Glukose ho.DPA_me.Glukose hormonomics DPA
## group2 to C.H
## ho.IAA_me.Glukose metabolomics Glukose -0.08555827
## ho.oxIAA_me.Glukose metabolomics Glukose 0.51707118
## ho.IAA.Asp_me.Glukose metabolomics Glukose -0.22176715
## ho.ABA_me.Glukose metabolomics Glukose 0.80887195
## ho.PA_me.Glukose metabolomics Glukose 0.86161713
## ho.DPA_me.Glukose metabolomics Glukose 0.73483118tail(e)## edge group1 from group2 to
## me.Val_qP.M0ZJG3 me.Val_qP.M0ZJG3 metabolomics Val qPCR M0ZJG3
## me.Met_qP.M0ZJG3 me.Met_qP.M0ZJG3 metabolomics Met qPCR M0ZJG3
## me.Ile_qP.M0ZJG3 me.Ile_qP.M0ZJG3 metabolomics Ile qPCR M0ZJG3
## me.Lys_qP.M0ZJG3 me.Lys_qP.M0ZJG3 metabolomics Lys qPCR M0ZJG3
## me.Leu_qP.M0ZJG3 me.Leu_qP.M0ZJG3 metabolomics Leu qPCR M0ZJG3
## me.Phe_qP.M0ZJG3 me.Phe_qP.M0ZJG3 metabolomics Phe qPCR M0ZJG3
## C.H
## me.Val_qP.M0ZJG3 0.9022587
## me.Met_qP.M0ZJG3 -0.3590730
## me.Ile_qP.M0ZJG3 0.8768180
## me.Lys_qP.M0ZJG3 0.8442707
## me.Leu_qP.M0ZJG3 0.5974332
## me.Phe_qP.M0ZJG3 0.4424355dim(e)## [1] 714 6# treatment 1
e1 <- extractEdges2(N1)
colnames(e1)[ncol(e1)] <- .treat[1]
head(e1)## edge group1 from
## ho.IAA_me.Glukose ho.IAA_me.Glukose hormonomics IAA
## ho.oxIAA_me.Glukose ho.oxIAA_me.Glukose hormonomics oxIAA
## ho.IAA.Asp_me.Glukose ho.IAA.Asp_me.Glukose hormonomics IAA.Asp
## ho.ABA_me.Glukose ho.ABA_me.Glukose hormonomics ABA
## ho.PA_me.Glukose ho.PA_me.Glukose hormonomics PA
## ho.DPA_me.Glukose ho.DPA_me.Glukose hormonomics DPA
## group2 to C
## ho.IAA_me.Glukose metabolomics Glukose -0.009706936
## ho.oxIAA_me.Glukose metabolomics Glukose 0.407917620
## ho.IAA.Asp_me.Glukose metabolomics Glukose -0.268229896
## ho.ABA_me.Glukose metabolomics Glukose 0.740096225
## ho.PA_me.Glukose metabolomics Glukose 0.896262542
## ho.DPA_me.Glukose metabolomics Glukose 0.829405769dim(e1)## [1] 662 6e <- merge(e,e1, sort=FALSE, all=TRUE)
head(e)## edge group1 from group2 to
## 1 ho.IAA_me.Glukose hormonomics IAA metabolomics Glukose
## 2 ho.oxIAA_me.Glukose hormonomics oxIAA metabolomics Glukose
## 3 ho.IAA.Asp_me.Glukose hormonomics IAA.Asp metabolomics Glukose
## 4 ho.ABA_me.Glukose hormonomics ABA metabolomics Glukose
## 5 ho.PA_me.Glukose hormonomics PA metabolomics Glukose
## 6 ho.DPA_me.Glukose hormonomics DPA metabolomics Glukose
## C.H C
## 1 -0.08555827 -0.009706936
## 2 0.51707118 0.407917620
## 3 -0.22176715 -0.268229896
## 4 0.80887195 0.740096225
## 5 0.86161713 0.896262542
## 6 0.73483118 0.829405769tail(e)## edge group1 from group2 to C.H
## 709 me.His_qP.HSP70 metabolomics His qPCR HSP70 0.8217718
## 710 me.Lys_qP.ACO2 metabolomics Lys qPCR ACO2 0.6368664
## 711 ho.oxIAA_me.Lys hormonomics oxIAA metabolomics Lys 0.2842323
## 712 me.Lys_qP.RbohA metabolomics Lys qPCR RbohA 0.5059654
## 713 me.Lys_qP.SnRK2 metabolomics Lys qPCR SnRK2 0.8284052
## 714 me.Lys_qP.P5CS metabolomics Lys qPCR P5CS 0.8220844
## C
## 709 NA
## 710 NA
## 711 NA
## 712 NA
## 713 NA
## 714 NA# treatment 2
.treat[2]## [1] "H"e2 <- extractEdges2(N2)
colnames(e2)[ncol(e2)] <- .treat[2]
head(e2)## edge group1 from
## ho.IAA_me.Glukose ho.IAA_me.Glukose hormonomics IAA
## ho.oxIAA_me.Glukose ho.oxIAA_me.Glukose hormonomics oxIAA
## ho.IAA.Asp_me.Glukose ho.IAA.Asp_me.Glukose hormonomics IAA.Asp
## ho.ABA_me.Glukose ho.ABA_me.Glukose hormonomics ABA
## ho.PA_me.Glukose ho.PA_me.Glukose hormonomics PA
## ho.DPA_me.Glukose ho.DPA_me.Glukose hormonomics DPA
## group2 to H
## ho.IAA_me.Glukose metabolomics Glukose -0.03135486
## ho.oxIAA_me.Glukose metabolomics Glukose 0.34316156
## ho.IAA.Asp_me.Glukose metabolomics Glukose -0.24733014
## ho.ABA_me.Glukose metabolomics Glukose 0.87529961
## ho.PA_me.Glukose metabolomics Glukose 0.88467893
## ho.DPA_me.Glukose metabolomics Glukose 0.89945125dim(e2)## [1] 688 6e <- merge(e,e2, sort=FALSE, all=TRUE)
head(e)## edge group1 from group2 to
## 1 ho.IAA_me.Glukose hormonomics IAA metabolomics Glukose
## 2 ho.oxIAA_me.Glukose hormonomics oxIAA metabolomics Glukose
## 3 ho.IAA.Asp_me.Glukose hormonomics IAA.Asp metabolomics Glukose
## 4 ho.ABA_me.Glukose hormonomics ABA metabolomics Glukose
## 5 ho.PA_me.Glukose hormonomics PA metabolomics Glukose
## 6 ho.DPA_me.Glukose hormonomics DPA metabolomics Glukose
## C.H C H
## 1 -0.08555827 -0.009706936 -0.03135486
## 2 0.51707118 0.407917620 0.34316156
## 3 -0.22176715 -0.268229896 -0.24733014
## 4 0.80887195 0.740096225 0.87529961
## 5 0.86161713 0.896262542 0.88467893
## 6 0.73483118 0.829405769 0.89945125tail(e)## edge group1 from group2 to C.H C H
## 735 me.Thr_qP.HSP70 metabolomics Thr qPCR HSP70 NA NA 0.2569343
## 736 me.Thr_qP.CAT1 metabolomics Thr qPCR CAT1 NA NA 0.3463942
## 737 me.Thr_qP.ERF1 metabolomics Thr qPCR ERF1 NA NA 0.1878546
## 738 me.Thr_qP.P5CS metabolomics Thr qPCR P5CS NA NA 0.2701877
## 739 me.Thr_qP.CO metabolomics Thr qPCR CO NA NA 0.3305915
## 740 me.Thr_qP.M0ZJG3 metabolomics Thr qPCR M0ZJG3 NA NA 0.3471071#write("Mid diablo 5 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!", "bla.log", append=TRUE)Compose file name and necessary information for network export file
file <- paste0("network-",paste(.treat, collapse="_"),"-",paste(datasets, collapse="_"),".txt")
label0 <- paste(paste(.treat, collapse=", "),"|",paste(datasets, collapse=", "),"; cutoff =",0)
title <- label0
sets <- 1:length(DATA)
suffix <- paste0(substr(names(DATA),1,2)[sets[-1]],collapse="-")write("Mid diablo 6 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!", "bla.log", append=TRUE)
write(file.path(.oroot,suffix,file), "bla.log", append=TRUE)Export edges table
if(!dir.exists(file.path(.oroot))) dir.create(file.path(.oroot))
if(!dir.exists(file.path(.oroot,suffix))) dir.create(file.path(.oroot,suffix))
#
my.write.table(e, file = file.path(.oroot,suffix,file)
, label=label0, na="0")## Warning in write.table(x, file = file, col.names = col.names, sep =
## sep, : appending column names to file## Object: e \\
## Label: C, H | hormonomics, metabolomics, qPCR ; cutoff = 0 \\
## File :\\
## \href{run:D:\\DEJAVNOSTI\\OMIKE\\pISA-projects\\multiOmics-integration\\_p_Omics\\_I_Omics\\_S_multiOmics\\_A_multiOmics-integration-R\\output\\ho-me-qP\\network-C_H-hormonomics_metabolomics_qPCR.txt}{../../../output/ho-me-qP/network-C_H-hormonomics_metabolomics_qPCR.txt}Table with edges for networks based on combined treatments (C, H) and single treatments (C) and (H) is exported as a text file. This table can be used for inspection and filtering out edges based on selected cutoff. Missing edges are labeled as weight 0. This enables numeric filtration in Excel.
write("End diablo 7 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!", "bla.log", append=TRUE)write("From 035-DIABLO !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!", "bla.log", append=TRUE)%%# Metadata files %% %%## Project %% %%{r pmeta, echo=FALSE,results='asis'} %%.pmeta %% %% %%## Investigation %% %%{r imeta, echo=FALSE,results='asis'} %%.imeta %% %% %%## Study %% %%{r smeta, echo=FALSE,results='asis'} %%.smeta %% %% %%## Assay %% %%{r ameta, echo=FALSE,results='asis'} %%.ameta %%
7 SessionInfo
Windows 10 x64 (build 19045)
R version 4.0.2 (2020-06-22) Platform: x86_64-w64-mingw32/x64 (64-bit) Running under: Windows 10 x64 (build 19045)
Matrix products: default
locale: [1] LC_COLLATE=Slovenian_Slovenia.1250 [2] LC_CTYPE=Slovenian_Slovenia.1250
[3] LC_MONETARY=Slovenian_Slovenia.1250 [4] LC_NUMERIC=C
[5] LC_TIME=Slovenian_Slovenia.1250
system code page: 1252
attached base packages: [1] grid stats graphics utils datasets grDevices [7] methods base
other attached packages: [1] pheatmap_1.0.12 ComplexHeatmap_2.6.2 igraph_1.2.6
[4] mixOmics_6.14.0 ggplot2_3.3.5 lattice_0.22-5
[7] MASS_7.3-60.0.1 pisar_0.1.0.9000 knitr_1.43
[10] rmarkdown_2.21
[4] tidyr_1.1.2 corpcor_1.6.9 png_0.1-7
[7] digest_0.6.35 RSpectra_0.16-0 R6_2.5.1
[10] plyr_1.8.6 stats4_4.0.2 ellipse_0.4.2
[13] evaluate_0.21 highr_0.8 pillar_1.4.7
[16] GlobalOptions_0.1.2 rlang_1.1.1 jquerylib_0.1.4
[19] S4Vectors_0.28.1 GetoptLong_1.0.5 Matrix_1.6-5
[22] labeling_0.4.2 rARPACK_0.11-0 stringr_1.4.0
[25] munsell_0.5.0 compiler_4.0.2 xfun_0.39
[28] pkgconfig_2.0.3 BiocGenerics_0.36.0 shape_1.4.6
[31] htmltools_0.5.2 tidyselect_1.1.0 tibble_3.0.4
[34] gridExtra_2.3 IRanges_2.24.1 matrixStats_1.2.0
[37] crayon_1.3.4 dplyr_1.0.2 withr_3.0.0
[40] jsonlite_1.8.8 gtable_0.3.0 lifecycle_0.2.0
[43] magrittr_2.0.1 scales_1.1.1 cli_3.2.0
[46] stringi_1.5.3 farver_2.0.3 reshape2_1.4.4
[49] bslib_0.3.1 ellipsis_0.3.1 generics_0.1.0
[52] vctrs_0.4.1 rjson_0.2.20 RColorBrewer_1.1-2 [55] tools_4.0.2 Cairo_1.5-15 glue_1.4.2
[58] purrr_0.3.4 parallel_4.0.2 fastmap_1.1.1
[61] yaml_2.2.1 clue_0.3-60 colorspace_1.4-1
[64] cluster_2.1.0 sass_0.4.0
%%prefix <- “”, sdir[i],suffix, “
Links to project on GitHub:
”) %%ind <- grep(“^_“,sdir) %%rind <- ind[1]-1 %% %%prefix <-”“, sdir[rind],suffix,” “,sep=”” %% ) %%for(i in 1:length(ind)) cat(” * “, %% file.path(prefix, sdir[rind],”tree/main”, paste(sdir[ind[1:i]], collapse=“/”)), “’>”, sdir[ind[i]],suffix, “”,sep=“” %% ) %% %%cat(“
“) %% %%# Link 1 %%cat(”
Main file :", mainFilePath, "
“) %%# %%mainFile <- strsplit(mainFile,‘.’,fixed=TRUE)[[1]][1] %%projectName <- rev((strsplit(dirname(getwd()),”/“))[[1]])[1] %%#cat(”Project file: [\href{run:“, %%#file.path(dirname(getwd()),paste0(projectName,”.prj”)), %%#“}{link}]
“,sep=”“) %%#